The Contribution of Epigenetics to Evolutionary Adaptation in Zingiber kawagoii Hayata (Zingiberaceae) Endemic to Taiwan

We epigenotyped 211 individuals from 17 Zingiber kawagoii populations using methylation-sensitive amplification polymorphism (MSAP) and investigated the associations of methylated (mMSAP) and unmethylated (uMSAP) loci with 16 environmental variables. Data regarding genetic variation based on amplified fragment length polymorphism (AFLP) were obtained from an earlier study. We found a significant positive correlation between genetic and epigenetic variation. Significantly higher mean mMSAP and uMSAP uHE (unbiased expected heterozygosity: 0.223 and 0.131, respectively, p < 0.001) per locus than that estimated based on AFLP (uHE = 0.104) were found. Genome scans detected 10 mMSAP and 9 uMSAP FST outliers associated with various environmental variables. A significant linear fit for 11 and 12 environmental variables with outlier mMSAP and uMSAP ordination, respectively, generated using full model redundancy analysis (RDA) was found. When conditioned on geography, partial RDA revealed that five and six environmental variables, respectively, were the most important variables influencing outlier mMSAP and uMSAP variation. We found higher genetic (average FST = 0.298) than epigenetic (mMSAP and uMSAP average FST = 0.044 and 0.106, respectively) differentiation and higher genetic isolation-by-distance (IBD) than epigenetic IBD. Strong epigenetic isolation-by-environment (IBE) was found, particularly based on the outlier data, controlling either for geography (mMSAP and uMSAP βE = 0.128 and 0.132, respectively, p = 0.001) or for genetic structure (mMSAP and uMSAP βE = 0.105 and 0.136, respectively, p = 0.001). Our results suggest that epigenetic variants can be substrates for natural selection linked to environmental variables and complement genetic changes in the adaptive evolution of Z. kawagoii populations.


Introduction
Adaptation to changing environments is a fundamental process for the survival of populations and species, especially during fast-paced environmental changes [1][2][3][4]. Epigenetics, which can occur without alterations in DNA sequences, is an important mechanism influencing population processes [5]. Over the ecological and evolutionary time scales of range redistribution, mutation and recombination would rarely provide sufficient sources of variation [6]. Epigenetic modifications have been suggested to respond to the environment before genetic changes begin to accumulate [7]. Phenotypic plasticity may arise as a result of epigenetic switches to deal with fluctuating environments [8], which may buy time for populations in the initial stages of adaptation [9,10]. Epigenetic modifications can generate heritable phenotypic variation in relation to adaptive evolution [11,12]. Population adaptive divergence in association with environmental gradients cannot be explained solely by DNA sequence variation [6]. Stable heritable epialleles may have significant evolutionary roles and are ecologically relevant in natural populations [13][14][15][16]. Epigenetic variation can be classified into (1) obligatory dependence of epigenetic variation on genetic variation, (2) facilitated epigenetic variation represents partial independence of epigenetic variation from genetic variation, and (3) pure epigenetic variation characterizes complete independence of epigenetic variation from genetic variation [17]. In natural plant populations, epigenetic variation can influence evolutionary processes in ways that are not related to sequence variation when epigenetic variation has arisen independently of genetic variation [6,17]. However, epigenetic variation may be the direct or indirect consequence of upstream genetic changes [6,[16][17][18]. Epigenetic variation, provoked by DNA methylation and demethylation, is an additional system compensating for adaptive genetic divergence [7,9,10]. Gene flow between populations may be reduced due to isolation-by-distance (IBD) [19]. IBD is the process by which geographically restricted gene flow generates a genetic structure indicating a positive correlation between genetic differentiation and geographic distance. The spatial genetic and epigenetic structure may differ, causing the difference between genetic and epigenetic IBD [20]. Additionally, ecological factors in divergent environments may lead to the selection-driven divergence that decreases gene flow between populations, creating a pattern of isolation-by-environment (IBE) [21,22].
In an earlier study, using data on amplified fragment length polymorphism (AFLP), a latitudinal pattern of environmentally dependent adaptive genetic divergence was found to be highly correlated with the annual temperature range in Z. kawagoii [23]. It is probable that the latitudinal northerly range expansion of Z. kawagoii [24] may have been related to the latitudinal pattern of adaptive divergence in Z. kawagoii [23]. Additionally, the lack of genetic variation in Z. kawagoii [23] may result in genetic as well as migration load [25,26]. Studies have demonstrated local adaptation linked to epigenetic variation and closely associated with environmental gradients in a variety of natural systems [27][28][29][30][31][32][33][34]. Therefore, the investigation of epigenetic variation in association with specific environmental variables is important. Epigenetic variation may play a role in compensation for the lack of genetic variation [7,9,10]. Testing for the association between epigenetic variation and environmental factors within a species' distribution range is important to identify environmental variables that may act as ecological drivers shaping the spatial epigenetic structuring of natural plant populations [6,16,27,[30][31][32][33][34]. The association between epigenetic variation and the environment may be related to alteration in the methylation status of different genes, resulting in local adaptation related to fitness-related traits [28,[34][35][36].
Epigenetic variation can be quantified with methylation-sensitive amplification polymorphism (MSAP) [37] that reflects a modification in cytosine methylation states [6,27,36]. In the present study, we surveyed the epigenetic variation in 211 individuals from 17 populations using MSAP to test the role that the environment plays in shaping adaptive epigenetic variation in Z. kawagoii. Epigenetic variation in association with environmental gradients can be important for plant adaptations to various environments apart from adaptive genetic divergence [5][6][7][9][10][11][12][13][14][15][16]. In addition to previous studies [13][14][15][16][27][28][29][30][31][32][33][34], the present study can provide additional empirical evidence for adaptive epigenetic divergence in natural plant populations. We aimed to answer questions related to spatial epigenetic structuring and environmentally dependent adaptive epigenetic divergence in Z. kawagoii. First, we explored the inter-population correlation between epigenetic distance and genetic distance using the Mantel test. Second, we intended to understand the level of population epigenetic diversity and epigenetic structure relative to the level of genetic diversity and genetic structure. Third, we assessed the correlations between all MSAP loci and environmental variables using a latent factor mixed model (LFMM) [38] and Samβada [39]. Fourth, F ST outliers were detected using genome scan methods including BAYESCAN [40] and DFDIST [41]. Fifth, F ST outliers identified using both BAYESCAN and DFDIST were further examined for their correlations with environmental variables using a Bayesian logistic regression approach [42]. Sixth, we assessed the linear fit of environmental variables with the ordination axes derived from a redundancy analysis (RDA) on the outlier epigenetic variation [43]. Subsequently, the environmental variables that showed a significant linear fit to the RDA ordinations were used in a partial RDA (pRDA) to examine the correlation between the environmental variables and pRDA axes conditioned on geographic effect. Lastly, we tested whether epigenetic IBE is present when controlling for either the geographic or genetic effect. The main objectives of the present study were to (1) assess the inter-population relationship between genetic and epigenetic variation in natural populations of Z. kawagoii, (2) understand how environmental variables influence population evolutionary epigenetic divergence and local adaptation, and (3) investigate if there is an environmentally dependent epigenetic divergence when controlling for the presence of any genetic structure.

Sample Collection and Epigenotyping
This study used the same 17 populations of Z. kawagoii examined in an earlier investigation [23] (Figure 1). Genomic DNA was extracted from silica gel dried leaf samples using a cetyltrimethyl ammonium bromide (CTAB) procedure [44], and ethanol-precipitated DNA was dissolved in 200 µL TE buffer (pH 8.0). The 211 plants used in this study included all but one individual from the TRK population used in the earlier AFLP investigation [23] due to a technical problem. A NanoDrop spectrophotometer (NanoDrop Technology, Wilmington, DE, USA) was used to quantify total genomic DNA. In brief, total genomic DNA (200 ng) was digested with the methylation-sensitive enzymes HpaII (1 U) and MspI (1 U) as frequent cutters separately with rare cutter EcoRI (1 U). Restriction digestion was performed in a total 10 µL reaction volume with 10X CutSmart buffer (New England Biolabs, Ipswich, MA, USA) for 1.5 h at 37 • C and then deactivated at 65 • C for 15 min. Restricted DNA products were ligated to MSAP adaptors (5 µM) with 5 U T4 DNA ligase (Thermo Scientific, Vilnius, Lithuania) and 5X ligation buffer (Thermo Scientific) at 22 • C for 1 h in a 10 µL ligation reaction mixture. using a Bayesian logistic regression approach [42]. Sixth, we assessed the linear fit of environmental variables with the ordination axes derived from a redundancy analysis (RDA) on the outlier epigenetic variation [43]. Subsequently, the environmental variables that showed a significant linear fit to the RDA ordinations were used in a partial RDA (pRDA) to examine the correlation between the environmental variables and pRDA axes conditioned on geographic effect. Lastly, we tested whether epigenetic IBE is present when controlling for either the geographic or genetic effect. The main objectives of the present study were to (1) assess the inter-population relationship between genetic and epigenetic variation in natural populations of Z. kawagoii, (2) understand how environmental variables influence population evolutionary epigenetic divergence and local adaptation, and (3) investigate if there is an environmentally dependent epigenetic divergence when controlling for the presence of any genetic structure.

Sample Collection and Epigenotyping
This study used the same 17 populations of Z. kawagoii examined in an earlier investigation [23] (Figure 1). Genomic DNA was extracted from silica gel dried leaf samples using a cetyltrimethyl ammonium bromide (CTAB) procedure [44], and ethanol-precipitated DNA was dissolved in 200 µL TE buffer (pH 8.0). The 211 plants used in this study included all but one individual from the TRK population used in the earlier AFLP investigation [23] due to a technical problem. A NanoDrop spectrophotometer (NanoDrop Technology, Wilmington, DE, USA) was used to quantify total genomic DNA. In brief, total genomic DNA (200 ng) was digested with the methylation-sensitive enzymes HpaII (1 U) and MspI (1 U) as frequent cutters separately with rare cutter EcoRI (1 U). Restriction digestion was performed in a total 10 µL reaction volume with 10X CutSmart buffer (New England Biolabs, Ipswich, MA, USA) for 1.5 h at 37 °C and then deactivated at 65 °C for 15 min. Restricted DNA products were ligated to MSAP adaptors (5 µM) with 5 U T4 DNA ligase (Thermo Scientific, Vilnius, Lithuania) and 5X ligation buffer (Thermo Scientific) at 22 °C for 1 h in a 10 µL ligation reaction mixture.   Table 1 for abbreviations of the population names. Pre-selective amplification was performed using a polymerase chain reaction (PCR). Adaptor ligated product (1:9 dilution with ddH 2 O) was mixed with 12 µM EcoRI primer (E00, Table S1), 12 µM HpaII-MspI primer (HM00, Table S1), 2.5 mM dNTPs, 1 U Taq DNA polymerase (Zymeset Biotech, Taipei, Taiwan), and 10X PCR buffer (Zymeset) in a 20 µL total volume. Pre-selective amplification was performed with an initial holding at 72 • C for 2 min and pre-denaturation at 94 • C for 3 min, followed by 25 cycles of 30 s at 94 • C, 30 s at 56 • C, and 1 min at 72 • C, with a final 5 min holding at 72 • C. Nine primer pair combinations with additional selective bases added at the ends of E00 and HM00 were used for selective amplification (Table S1). Fluorescent-dye-labeled 10 µM EcoRI selective primer was mixed with 10 µM HpaII-MspI primer, 2.5 mM dNTPs, 2 U Taq DNA polymerase (Zymeset), 10X PCR buffer (Zymeset), and diluted pre-selective amplified products (1:19 dilution with ddH 2 O) in a 20 µL total volume. Subsequently, selective PCR was performed with an initial holding at 94 • C for 3 min, followed by 13 cycles of 30 s at 94 • C, 30 s at 65-56 • C (decreasing the temperature by 0.7 • C each cycle), 1 min at 72 • C, and then 23 cycles of 30 s at 94 • C, 30 s at 56 • C, and 1 min at 72 • C, with a final 5 min holding at 72 • C. Selective amplification products were electrophoresed with an ABI 3730XL DNA analyzer (Applied Biosystem, Foster City, CA, USA).
We scored amplification fragments in the range of 150-500 bp, setting a fluorescent threshold of 150 units using Peak Scanner v.1.0 (Applied Biosystem). Low peak fragments and those scored higher than 99% or less than 1% of individuals were removed. Additionally, fragments within one nucleotide in a ± 0.8 base pair window were recognized as the same fragment. To determine the DNA methylation status of every locus, the "mixed scoring 1" transformation scheme in the R function MSAP-calc [45] in the R environment [46] was applied to obtain mMSAP (methylated) and uMSAP (unmethylated) datasets. Three randomly chosen samples in each population for each primer combination were used to assess the epigenotyping error rate. The error rate for EcoRI-MspI (eMspI), EcoRI-HpaII (eHpaII), and a combined error rate (eMspI + eHpaII − 2e MspIeHpaII ) for each primer combination was calculated [27]. We removed loci with an error rate per locus greater than 5% [47]. The mean error rate was 2.15 and 2.12%, respectively, for eMspI and eHpaII, with a combined error rate of 4.18% (Table S1).

Epigenetic Diversity, Differentiation, and Clustering
The software AFLP-SURV v.1.0 [48] was used to estimate population unbiased expected heterozygosity (uH E ) [49] and the percentage of polymorphic loci (PPL) based on allele frequencies using the settings of the Hardy-Weinberg equilibrium and a non-uniform prior distribution. Per locus uH E was estimated using ARLEQUIN v.6.0 [50]. A linear mixed effect model (LMM) was used to estimate the difference in mean uH E per locus between three types of markers (AFLP, mMSAP, and uMSAP). AFLP data were obtained from an earlier study [23]. In LMM, the marker type and population were used as a fixed factor and a random factor, respectively, and analyzed using the lmer function in the R package lme4 [51]. Significance was tested using the Anova function in the R package car based on type II Wald χ 2 statistic [52]. Tukey's post hoc pairwise comparisons of mean uH E per locus between marker types were assessed using the lsmeans function in the R package emmeans [53]. Additionally, a paired t-test was used to assess the differences in mean uH E between the three marker types at the population level. Population genetic (AFLP) and epigenetic (mMSAP and uMSAP) distance matrices were calculated using Nei's genetic distance [54] using the nei.dist function in the R package poppr [55]. Mantel correlations for the inter-population relationship between genetic and epigenetic variation, based on these distance matrices, were assessed using the mantel function in the R package vegan [56]. A partial Mantel test, performed using the mantel.partial function in the R package vegan, was also used to assess the inter-population relationship between genetic and epigenetic variation controlling for the geographic effect. The geographic effect was computed using the coordinates of the sample sites.
Analysis of molecular variance (AMOVA) was used to estimate the level of genetic differentiation between populations (Φ ST ) using the poppr.amova function in the R package poppr. Significance was tested using the randtest function in the R package ade4 [57] with 9999 permutations. Pairwise population F ST was computed using ARLEQUIN, and significance was tested with 10,000 permutations. Pairwise population F ST values were used to calculate the level of divergence for each population from the remaining populations as the mean value of pairwise F ST for each population against the rest of the populations (denoted as the population mean F ST ). Epigenetic homogeneous groups of individuals were assessed using discriminant analysis of principal components (DAPC) [58]. The find.clusters and dapc functions in the R package adegenet [59] were used in DAPC analysis setting K = 1-10. The Bayesian information criterion (BIC) in DAPC was estimated to determine the optimal number of clusters.

Test for F ST Outliers
To detect signatures in selection on MSAP loci, F ST outliers were identified using BAYESCAN and DFDIST. BAYESCAN v.2.1 [40] was used to estimate the ratio of posterior probabilities of selection over neutrality (the posterior odds (PO), log 10 (PO)). Two hundred pilot runs of 50,000 iterations followed by a sample size of 50,000 with a thinning interval of 20 among 10 6 iterations were performed in BAYESCAN. A logarithmic scale of log 10 (PO) > 1 was used as strong evidence (posterior probability > 0.91) for selection over neutrality for a locus under directional selection [60]. DFDIST estimates a distribution of observed F ST versus uH E , and loci under selection were identified by comparing them to a simulated neutral distribution. In DFDIST, we set critical frequency = 0.99, Zhivotovsky parameter = 0.25, trimmed mean F ST = 0.3 (excluding 30% of highest and 30% of lowest F ST values), smoothing proportion = 0.04, 500,000 resamplings, and critical p = 0.05. MSAP loci with observed F ST against uH E falling above the 95% confidence level of the simulated distribution were recognized as F ST outliers under directional selection.

Test for Associations between Epigenetic Loci and Environmental Variables
The associations between all epigenetic loci and environmental variables were estimated using LFMM and Samβada. In LFMM, we considered population epigenetic structure as a random factor using a latent factor of 1 and 4, respectively, according to the DAPC results (see Results) for mMSAP and uMSAP. Matrices of mMSAP and uMSAP variation were used as fixed factors. Ten LFMM runs with 10,000 iterations of the Gibbs sampling algorithm and a burn-in period of 5000 cycles were performed for each environmental predictor. Z-scores for each environmental predictor were obtained by combining the results of ten independent LFMM runs, and p-values were adjusted using the genomic inflation factor [38]. Moreover, a 1% false discovery rate (FDR) was used to adjust the p-value using the qvalue function in the R package qvalue [61]. A multiple univariate logistic regression approach in Samβada was used to assess the correlations between the allele frequencies of mMSAP and uMSAP loci and the values of the environmental variables. Both Wald and G scores [39] were used, applying a 1% FDR for the p-value adjustment, to assess the fit of the model with environmental variables against the null model without environmental variables.
Moreover, a Bayesian logistic regression analysis, implemented with the stan_glm function in the R package rstanarm [42], was used to further justify the associations between environmental variables and the potential F ST outliers that were identified with both BAYESCAN and DFDIST. The weakly informative priors following a student's t distribution with mean zero and seven degrees of freedom were used in stan_glm, and the scale of the prior distribution was 10 for the intercept and 2.5 for the predictors. We ran all the stan_glm models using 4 chains, each containing 2000 warm-up and 2000 sampling steps, and a 95% credible interval was determined using the posterior_interval function in the R package rstanarm. In the stan_glm analysis, we obtained values for the effective sample size and convergence diagnostic statistic representing good priors applied and stable estimates obtained for each predictor.

Linear Relationships between Environmental Variables and the Ordination Axes of the Redundancy Analysis
A multivariate approach in an RDA analysis [43] was used to estimate the extent to which outlier variation was explained by the environmental variables. We explored the relationships between epigenetic variation and environmental drivers by fitting variables onto ordinations using the envfit function in the R vegan package. In the envfit analysis, a full RDA model was used to estimate the independent effect of the environment (16 environmental variables) fitting to the amount of outlier variation with 999 permutations and represented by the squared correlation coefficients (R 2 ). All p-values were adjusted for multiple comparisons using 5% FDR. Environmental variables with significant fit to the outlier epigenetic variation were then used in a pRDA analysis. pRDA was analyzed in order to assess the correlations between environmental variables and the first two axes conditioned on the geographic effect. Significance of the pRDA ordination axes was assessed using the anova.cca function in the R package vegan with 999 permutations. The arrows pointed in the direction of maximum variation in the value of environmental variables, and the degree to which the variables correlated with pRDA axes was represented by the length of the arrows.

Epigenetic Isolation-by-Environment Controlling for Geographic or Genetic Effects
The Mantel function in the R package vegan (999 permutations) was used to test IBD based on the genetic and epigenetic distance matrices against the geographic distance Plants 2023, 12, 1558 7 of 18 matrix calculated using the dist function in R. To test IBE, Mantel and partial Mantel tests (controlling for the geographic effect) were used to assess the relationship between the epigenetic and environmental distance matrices, respectively, using the mantel and mantel.partial functions implemented in the R package vegan. IBD and IBE were also performed using the MMRR (multiple matrix regression with randomization) function in R [21]. In addition, IBE was also tested for epigenetic distance matrices against the environmental distance matrix controlling for genetic effect (AFLP) using the partial Mantel test and MMRR. This was performed to test if environmental variables influence epigenetic variation independent of the genetic effect (isolation-by-genetic structure, IBG). In MMRR, regression coefficients for IBD (β D ), IBG (β G ), and IBE (β E ) were obtained, and the significance was determined after 999 permutations.

Epigenetic Diversity and Structure
Overall, 9 primer pairs resolved a total of 481 unambiguous bands ranging from 22 to 107 (Table S1), and the number of loci estimated using the R MSAP-calc script [45] were 424 and 354, respectively, for mMSAP and uMSAP. The average PPL was 58.8% in mMSAP and 33.6% in uMSAP. The average uH E was 0.160 in mMSAP and 0.119 in uMSAP ( Table 1). The LMM analysis showed a significant difference in mean uH E per locus among the three types of markers (AFLP, mMSAP, and uMSAP; χ 2 = 2002.8, p < 0.0001). Pairwise comparisons between marker types revealed that the mean uH E per locus was significantly smaller in AFLP (mean uH E per locus = 0.104) compared to mMSAP (mean uH E per locus = 0.223) and uMSAP (mean uH E per locus = 0.131) (Table S3). At the population level, genetic uH E was not significantly different from uMSAP epigenetic uH E (paired t-test: t 16 = 0.796, p = 0.438), but it was significantly lower compared to mMSAP epigenetic uH E (paired t-test, t 16 = −7.426, p < 0.0001). mMSAP uH E was significantly higher compared to uMSAP uH E (paired t-test, t 16 = 8.007, p < 0.0001).
The overall values of F ST based on all loci were 0.044 and 0.106, respectively, for mMSAP and uMSAP (Table S4). The AMOVA showed an overall population differentiation (Φ ST ) of 0.071 and 0.181, respectively, for mMSAP and uMSAP based on the total data ( Table 2). The average F ST values computed were 0.044 and 0.106, respectively, for the total mMSAP and uMSAP variation (Table S4). Using the total data, the lowest BIC was found at K = 2 and K = 6, respectively, for mMSAP and uMSAP in DAPC ( Figure S1). No distinct mMSAP population clustering was found (Figure 2a), but four uMSAP population clusters were revealed (Figure 2b). However, individuals in the same population may be grouped in different uMSAP clusters. Unlike AFLP [23], the annual temperature range and latitude showed no significant relationship with the mMSAP population mean

Inter-Population Relationship between Genetic and Epigenetic Variation
A pairwise population Nei's distance matrix was used in a Mantel test to investigate the inter-population correlation between genetic and epigenetic variation. We found significant inter-population correlations between genetic and epigenetic distance (AFLP vs. mMSAP: Mantel r = 0.484, p = 0.002; AFLP vs. uMSAP: Mantel r = 0.596, p = 0.001; and mMSAP vs. uMSAP: Mantel r = 0.832, p = 0.001), and the linear fitting lines are displayed in Figure 3. The relationships between the three marker types were also significant when the effect of geography was excluded using the partial Mantel test (AFLP vs. mMSAP: Mantel r = 0.388, p = 0.007; AFLP vs. uMSAP: Mantel r = 0.460, p = 0.004; and mMSAP vs. uMSAP: Mantel r = 0.814, p = 0.001). Table 2. Epigenetic differentiation between the 17 populations of Zingiber kawagoii based on the total and outlier mMSAP and uMSAP variation using an analysis of molecular variance (AMOVA).    nificant inter-population correlations between genetic and epigenetic distance (AFLP vs. mMSAP: Mantel r = 0.484, p = 0.002; AFLP vs. uMSAP: Mantel r = 0.596, p = 0.001; and mMSAP vs. uMSAP: Mantel r = 0.832, p = 0.001), and the linear fitting lines are displayed in Figure 3. The relationships between the three marker types were also significant when the effect of geography was excluded using the partial Mantel test (AFLP vs. mMSAP: Mantel r = 0.388, p = 0.007; AFLP vs. uMSAP: Mantel r = 0.460, p = 0.004; and mMSAP vs. uMSAP: Mantel r = 0.814, p = 0.001).

FST Outlier Identification and Association between Outliers and Environmental Variables
In the DFDIST analysis, 41 mMSAP and 10 uMSAP loci (9.67 and 2.82%, respectively) were identified as being FST outliers ( Figure S2). The BAYESCAN analysis identified 45 mMSAP and 15 uMSAP loci as being FST outliers, corresponding to 10.61 and 4.24% of all loci, respectively. Although we set log10(PO) > 1.0 as a criterion for strong evidence in the identification of outliers using BAYESCAN, most of the outliers identified with BAYESCAN had a log10(PO) > 2 (decisive evidence for selection corresponding to a probability > 0.99) [58], except two mMSAP outliers (mX06HM_2011 and mX10HM_4950) (Table S5). A total of 10 mMSAP and 9 uMSAP loci, respectively, were identified as outliers using both DFDIST and BAYESCAN (Table S5, Figure S2). No outlier mMSAP loci were correlated with BIO19 or CLO assessed using LFMM, Samβada, and rstanarm (Table S5). Most outlier loci exhibited a significant association with two or more environmental variables based on the three regression approaches (Table S5). The ΦST was 0.284 in mMSAP and 0.404 in uMSAP based on the outlier loci (Table 2).

Environmental Effect on Outlier Epigenetic Variation
Using the 16 environmental variables in the full RDA model, the statistical test supported the role of the environment in shaping the distribution of the outlier mMSAP and uMSAP epigenotypes (mMSAP: adjusted R 2 = 0.273, F = 5.917, p = 0.001; uMSAP: adjusted R 2 = 0.390, F = 9.414, p = 0.001). The full RDA model analysis with envfit suggested that outlier mMSAP and uMSAP variation, respectively, correlated significantly with 11 and 12 environmental variables (5% FDR adjusted p < 0.05, Table 3). These 11 and 12 environmental variables were used, respectively, in the pRDA model for outlier mMSAP and

F ST Outlier Identification and Association between Outliers and Environmental Variables
In the DFDIST analysis, 41 mMSAP and 10 uMSAP loci (9.67 and 2.82%, respectively) were identified as being F ST outliers ( Figure S2). The BAYESCAN analysis identified 45 mM-SAP and 15 uMSAP loci as being F ST outliers, corresponding to 10.61 and 4.24% of all loci, respectively. Although we set log 10 (PO) > 1.0 as a criterion for strong evidence in the identification of outliers using BAYESCAN, most of the outliers identified with BAYESCAN had a log 10 (PO) > 2 (decisive evidence for selection corresponding to a probability > 0.99) [58], except two mMSAP outliers (mX06HM_2011 and mX10HM_4950) (Table S5). A total of 10 mMSAP and 9 uMSAP loci, respectively, were identified as outliers using both DFDIST and BAYESCAN (Table S5, Figure S2). No outlier mMSAP loci were correlated with BIO19 or CLO assessed using LFMM, Samβada, and rstanarm (Table S5). Most outlier loci exhibited a significant association with two or more environmental variables based on the three regression approaches (Table S5). The Φ ST was 0.284 in mMSAP and 0.404 in uMSAP based on the outlier loci (Table 2).

Environmental Effect on Outlier Epigenetic Variation
Using the 16 environmental variables in the full RDA model, the statistical test supported the role of the environment in shaping the distribution of the outlier mMSAP and uMSAP epigenotypes (mMSAP: adjusted R 2 = 0.273, F = 5.917, p = 0.001; uMSAP: adjusted R 2 = 0.390, F = 9.414, p = 0.001). The full RDA model analysis with envfit suggested that outlier mMSAP and uMSAP variation, respectively, correlated significantly with 11 and 12 environmental variables (5% FDR adjusted p < 0.05, Table 3). These 11 and 12 environmental variables were used, respectively, in the pRDA model for outlier mMSAP and uMSAP conditioned on geography. The first two axes of the pRDA explained 68.81 and 73.75% of the outlier mMSAP and uMSAP variation, respectively (mMSAP: adjusted R 2 = 0.140, F = 4.427, p = 0.001; uMSAP: adjusted R 2 = 0.269, F = 8.224, p = 0.001) (Figure 4). Table 3. Environmental covariates that were significantly associated with outlier epigenetic variation based on the envfit analysis. R 2 , the amount of variation explained by each covariate in the model, was determined using envfit. Adjusted p, the p-value corrected for multiple comparisons, was determined using 5% FDR after 999 permutation tests. uMSAP conditioned on geography. The first two axes of the pRDA explained 68.81 and 73.75% of the outlier mMSAP and uMSAP variation, respectively (mMSAP: adjusted R 2 = 0.140, F = 4.427, p = 0.001; uMSAP: adjusted R 2 = 0.269, F = 8.224, p = 0.001) (Figure 4).  When conditioned on the geographic effect, the outlier mMSAP variation was highly correlated with slope, soil pH, and BIO12 on the first pRDA axis (Table 4 and Figure 4). When conditioned on the geographic effect, the outlier mMSAP variation was highly correlated with slope, soil pH, and BIO12 on the first pRDA axis (Table 4 and Figure 4). Aspect, NDVI, PET, and BIO12 were highly correlated with the outlier mMSAP variation on the second pRDA axis. On the first pRDA axis, the outlier uMSAP variation showed strong correlations with aspect, slope MI, soil pH, and BIO12. The second pRDA axis revealed high correlations between MI and NDVI and the outlier uMSAP variation. The annual temperature range had a relatively higher R 2 (mMSAP: R 2 = 0.308; uMSAP: R 2 = 0.129) than the other environmental variables in the envfit analysis (Table 3). However, the annual temperature range had low levels of correlation with the outlier mMSAP (RDA1 biplot score = 0.074, RDA2 biplot score = −0.010) and uMSAP (RDA1 biplot score = 0.019, RDA2 biplot score = −0.029) variation in both pRDA axes when controlling for the geographic effect (Table 4 and Figure 4).

Discussion
Theoretical work has shown that beneficial and heritable epigenetic variants may play important roles for populations to adapt quickly to local conditions independent of genetic changes [7]. Although empirical studies have shown epigenetics may act independently from underlying genetic variation [29,30,33,62,63], genotypes and epigenotypes may interact and together affect biological processes, for example, in Viola cazorlensis [27], Arabidopsis thaliana [28], and Taiwania cryptomerioides [64]. Our results indicate inter-population correlations between genetic and epigenetic variation, even with the exclusion of the geographic effect using the partial Mantel test (Figure 3). The underlying mechanism could be that the epigenetic variation interplay with genetic variation and epigenetic variation is at least in part a downstream, subsidiary effect of genetic changes [6,[16][17][18]. Epigenetic variation may have played a role in speeding up the local adaptation of Z. kawagoii. Relatively lower Z. kawagoii AFLP diversity was found compared to that of Zingiber species distributed in Brazil and India [23]. We found that mean mMSAP diversity is relatively higher compared to mean AFLP diversity (Table 1, cf. 23). At the population level, most populations of Z. kawagoii had a relatively higher level of epigenetic than genetic diversity, particularly when compared between AFLP and mMSAP ( Table 1, cf. 23). Additionally, mean uH E per locus was significantly smaller in AFLP compared to mMSAP and uMSAP (Table S3). These results suggest that epigenetic mutation may occur at a faster rate than genetic mutation [65,66]. Moreover, lower epigenetic than genetic differentiation (Table 2 and  Table S4; cf. 23) suggests larger epigenetic variation in shorter spatial distances might contribute considerably to population epigenetic diversity [27,33,34,67,68]. It is possible that epigenetic variation may provide a compensating source for the low level of genetic variation [7,9,10].
Neutral epigenetic divergence may occur via drift, and inter-population correlation between genetic and epigenetic variation may involve no functional link [5,69]. The higher genetic than epigenetic differentiation found in Z. kawagoii and other plants [30,63,64,70] ( Table 2 and Table S4; cf. 23) suggests that geographic distance was a better predictor of genetic than epigenetic differentiation. This is evidenced by the larger Mantel statistic for AFLP relative to mMSAP and uMSAP in this study (AFLP: Mantel r = 0.457; mMSAP: Mantel r = 0.126; uMSAP: Mantel r = 0.184). Nonetheless, significant epigenetic IBD was detected with both the Mantel test and MMRR based on the total and outlier data ( Table 5), suggesting that population epigenetic variation may be partly caused by random epigenetic drift [71]. Additionally, the stronger signal in genetic IBD than in epigenetic IBD suggests that the signal is biased by non-unidirectional changes in environmental gradients, which may impede epigenetic IBD [72,73].
Epigenetic changes can be triggered by biotic and abiotic stresses in natural environments and provide potential for novel heritable variation [74]. Environmental changes affect epigenetic variation in many organisms, which might help them adapt to rapid environmental changes [27][28][29][30]33]. Thus, the role of epigenetics in response to environmental stimulus has been an important issue in evolution [36]. With a significant inter-population correlation between genetic and epigenetic variation, strong epigenetic IBE can still be found when controlling for the effect of any genetic structure using both the partial Mantel test and MMRR (Table 5), particularly based on the outlier data. It is likely that part of the epigenetic variation is not coupled with genetic changes, and environments may act as causal factors in inducing epigenetic variation. Thus, epigenetic variation may play a role in the rapid evolution of individuals to various environments [64]. Epigenetic profiles may diverge between environments contributing to differentially induced effects [29,75]. Additionally, there was a significantly lower level of mean genetic uH E per locus compared to mean epigenetic uH E per locus (Table S3), reinforcing the role of epigenetic variation in Z. kawagoii adaptation. This may enable higher survival in a population encountering novel environments until the acquisition of beneficial genetic mutation.
The current study extends the earlier research [23] by examining the links between environmental variables and epigenetic variants. Unlike the earlier study [23], which showed that annual temperature range was the main driver for latitudinal adaptive divergence, this study found that the annual temperature range may play only a minor role in the adaptive epigenetic divergence, as suggested by low biplot scores in the first and second pRDA axes (Table 4 and Figure 4). This suggests low correlation between the annual temperature range and the outlier MSAP variation at larger spatial scales. Moreover, epigenetics in this study, similar to the earlier genetic study [23], indicated that combinations of different environmental drivers may act as selective pressures in invoking adaptive evolution. However, environmental variables may have differential effects on adaptive genetic and epigenetic variation in Z. kawagoii. The results of pRDA analyses in this study suggest outlier mMSAP variation is associated mainly with slope, soil pH, NDVI, PET, and BIO12, whereas outlier uMSAP variation is mainly associated with aspect, slope, MI, NDVI, soil pH, and BIO12 (Table 4 and Figure 4). However, the contributions of various environmental variables with high correlations with both pRDA axes (Table 4 and Figure 4) suggest that these variables may have acted in concert on various genes invoking adaptive epigenetic divergence.
Although it is unclear whether environmentally dependent epigenetic changes at distinct loci are the direct consequence of environments, local environments may play a  [36,72,76]. This may be advantageous if epigenetic changes that are stress responsive and are ecologically relevant. For example, individuals of L. racemosa from salt marsh populations display a shrub-like phenotype with lower levels of DNA methylation, whereas individuals of riverside populations with a tree-like phenotype had higher levels of DNA methylation [34]. In this mangrove species, genes may be activated via demethylation due to salinity stress in the salt marsh populations in contrast to hypermethylation of these genes in the riverside populations. Correlations between epigenetic changes and traits related to plant growth and development can be found, for example, leaf traits in Prunus mume [77] and Carpobrotus edulis [78], flower morphology in V. cazorlensis [27], and inhabiting different natural habitats in Phragmites australis and Lilium bosniacum [79,80]. Epigenetic changes are closely related to climatic factors such as temperature [31,33,81] and precipitation [31,33]. Epigenetic changes are also found to be closely associated with ecological factors such as soil factors [31,82], NDVI, and PET [31,33]. Topographic variables including slope and aspect are also found to be closely related to epigenetic changes [31,33]. Therefore, environmentally dependent epigenetic changes may be related to the growth and development of plants, thus contributing to ecological and evolutionary processes in natural populations.

Conclusions
Genetic adaptation to environmental changes is crucial for the conservation and survival of species. In addition, researchers in evolutionary biology are increasingly interested in epigenetic adaptations to different environments. Both genetic and epigenetic variations may be sources of variation that play important roles in adapting to environmental heterogeneity and, hence, are important for understanding how the environment shapes natural population diversity in a changing environment. In the present study, epigenetic variation provided no distinct regional substructuring, particularly based on the total mMSAP data. This study revealed low epigenetic population differentiation, indicating that larger epigenetic variation occurs mainly at shorter spatial scales. Epigenetic variation may have compensated for the lack of genetic variation and played an important role in the adaptive divergence and local adaptation in Z. kawagoii populations. We identified that outlier MSAP loci associated strongly with specific environmental variables that may have played important roles in the local adaptation and survival of Z. kawagoii populations despite significant IBD. The local adaptation involving epigenetic changes in Z. kawagoii is also evidenced by the finding of significant IBE based on the outlier epigenetic variation when controlling for the geographic or genetic effects. Nonetheless, this study suggests that selection may be involved in population divergence at epigenetic sites that are partly linked to genetic sites in Z. kawagoii. Our findings highlight that methylation changes can be substrates for natural selection associated with environmental variables and may have a long-term consequence on the adaptation and survival of Z. kawagoii populations. This study also emphasizes the importance of investigating population epigenetic variation in association with environmental gradients apart from adaptive genetic divergence.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/plants12071558/s1, Figure S1: Bayesian information criterion for the evaluation of clustering scenarios analyzed based on mMSAP and uMSAP. Figure. S2: Venn diagrams showing the number of outlier MSAP (mMSAP and uMSAP) loci identified using genome scans of natural populations of Zingiber kawagoii using BAYESCAN and DFDIST. Table S1: Primer combinations, number of markers, and error rate per locus calculated using the MSAP technique.  Table S3: Summary of Tukey's post hoc pairwise comparisons of the mean unbiased expected heterozygosity (uH E ) per locus in Zingiber kawagoii between marker types (AFLP, mMSAP, and uMSAP) using a linear mixed effect model. Table S4: Matrix of pairwise F ST (lower diagonal) and p-values (upper diagonal) based on the total mMSAP and uMSAP data of the 17 populations of Zingiber kawagoii estimated using ARLEQUIN. Table S5: F ST outliers identified using BAYESCAN and DFDIST that were strongly associated with environmental variables. BAYESCAN and DFDIST. Table S1: Primer combinations, number of markers, and error rate per locus calculated using the MSAP technique.  Table S3: Summary of Tukey's post hoc pairwise comparisons of the mean unbiased expected heterozygosity (uHE) per locus in Zingiber kawagoii between marker types (AFLP, mMSAP, and uMSAP) using a linear mixed effect model. Table S4: Matrix of pairwise FST (lower diagonal) and p-values (upper diagonal) based on the total mMSAP and uMSAP data of the 17 populations of Zingiber kawagoii estimated using ARLEQUIN. Table S5: FST outliers identified using BAYESCAN and DFDIST that were strongly associated with environmental variables.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.